##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## R code for various functions needed for the replication of the paper
## "Incumbency Advantages, Prime Minister Replacements and Government Formation"
## by Hanna Bäck, Marc Debus, and Michael Imre
##
## the original likelihood functions were written by Chiba et al. (2015) and adapted by Ecker and Meyer (2020) (see paper for full citations)
## the optimization in estimate.cld was further adapted by Michael Imre.
## hmftest.clogit was adapted from hmftest (from the mlogit package) by Michael Imre.




# log-likelihood functions for univariate models --------------------------

## Logit model

#logit.ll <- function(vec.betas, data){
#  loglik <- sum(-data$y.vec*log(1 + exp(-(data$x.mat%*%vec.betas))) - (1-data$y.vec)*log(1 + exp(data$x.mat%*%vec.betas)))
#  return(-loglik)
#}

logit.ll.prob <- function(vec.betas, data){
  ## Univariate logit estimator.
  ## It returns log-likelihood values given data and parameters.
  ##  pars = vector of parameters
  ##  data = list of 1 vector and 1 matrix.
  ##   (1) selection var (y.vec),
  ##   (2) selection covariates (x.mat)
  #prob <- exp(data$x.mat%*%vec.betas)/(1+exp(data$x.mat%*%vec.betas))
  prob <- 1/(1+exp(-data$x.mat%*%vec.betas))
  -sum((data$y.vec*log(prob))+((1-data$y.vec)*log(1-prob)))
  #loglik <- sum((data$y.vec*log(prob))+((1-data$y.vec)*log(1-prob)))
  #return(-loglik)
}

## Conditional logit model

clogit.ll <- function(pars, data, const.vec = 0){
	## Univariate c-logit estimator.
	## It returns log-likelihood values given data and parameters.
	##  pars = vector of parameters
	##  data = list of 2 vectors and 1 matrix.
	##   (1) selection var (y.sel),
	##   (2) grouping var (caseid)
	##   (3) duration covariates (x.sel)
	agg <- exp(data $ xmat %*% pars + const.vec) ## add constraints to the aggregator
	fac <- factor(data $ caseid)
	denom <- tapply(agg, fac, sum) ## only for y==1
	prob <- agg[data $ y==1]/denom
	prob[is.na(prob)] <- 1e-20
	sum(log(prob))
}

## Weibull duration model

weibull.ll <- function(pars, data){
	## Univariate Weibull
	## It returns log-likelihood values given data and parameters.
	##  pars = vector of parameters
	##  data = list of 2 vectors and 1 matrix.
	##   (1) duration var (y.dur),
	##   (2) right censoring indivator (t=0 means censoring)
	##   (3) duration covariates (x.dur)
	k <- ncol( data $ xmat )
	p <- exp( pars[k+1] )
	lambda <- exp( -data$xmat %*% pars[1:k] )
	logft <- log( p * lambda )+( (p-1) * log( data$y * lambda) ) -( data$y * lambda )^p
	logSt <-  -( data $ y * lambda )^p
	loglik <- (data$t==1)*logft + (data$t==0)*logSt
	sum(loglik)
}

## Weibull stratified duration model

weibull.ll.str <- function(pars, data){
  ## Univariate Weibull
  ## It returns log-likelihood values given data and parameters.
  ##  pars = vector of parameters
  ##  data = list of 2 vectors and 1 matrix.
  ##   (1) duration var (y.dur),
  ##   (2) right censoring indivator (t=0 means censoring)
  ##   (3) duration covariates (x.dur)
  k <- ncol( data $ xmat )
  p0 <- pars[k+1]
  p1 <- pars[k+1] %*% strata
  p <- exp( p0 )
  lambda <- exp( -data$xmat %*% pars[1:k] )
  logft <- log( p * lambda )+( (p-1) * log( data$y * lambda) ) -( data$y * lambda )^p
  logSt <-  -( data $ y * lambda )^p
  loglik <- (data$t==1)*logft + (data$t==0)*logSt
  sum(loglik)
}

## Log-logistic duration model

loglogistic.ll <- function(pars, data){
	## Univariate Log-Logistic
	## It returns log-likelihood values given data and parameters.
	##  pars = vector of parameters
	##  data = list of 2 vectors and 1 matrix.
	##   (1) duration var (y.dur),
	##   (2) right censoring indivator (t=0 means censoring)
	##   (3) duration covariates (x.dur)
	k <- ncol( data$xmat )
	p <- exp( pars[k+1] )
	lambda <- exp( - data $ xmat %*% pars[1:k] )
	logft <- log(lambda) + log(p) + (p-1) * log( data $ y * lambda) - 2 * log(1+ ( data $ y * lambda )^p)
	logSt <-  -log( 1 + ( data $ y * lambda )^p )
	loglik <- (data$t==1)*logft + (data$t==0)*logSt
	sum(loglik)
}

# log-likelihood functions for bivariate models ---------------------------

## Weibull duration with logit selection

log.weib.ll <- function(pars, data){
  ## Bivariate c-logit-Weibull estimator for FIML estimation.
  ## It returns log-likelihood values given data and parameters.
  ##	pars = vector of parameters
  ##	data = list of 4 vectors and 2 matrices.
  ##   (1) selection var (y.sel)        (2) duration var (y.dur),
  ##   (3) selection covariates (x.sel) (4) duration covariates (x.dur)
  ##   (5) grouping var (caseid)        (6) right censoring indivator (t=0 means censoring)
  ##==========================##
  ## First, extract parameters
  ##==========================##
  k.sel <- ncol(data $ x.sel) # n(ind v) for selection
  k.dur <- ncol(data $ x.dur) # n(ind v) for outcome
  b.sel <- pars[1:k.sel]
  b.dur <- pars[(k.sel+1):(k.sel+k.dur)]
  p <- exp(pars[(k.sel+k.dur+1)]) # log(p) -> p
  b.rho <- (exp(2*pars[length(pars)])-1)/(exp(2*pars[length(pars)])+1) # athrho -> rho
  b.rho <- min(b.rho, .9999, na.rm=T) # This is to avoid denom=0
  t <- data$t	# censoring indicator
  ##===============##
  ## logit part
  ##===============##
  #prob <- 1/(1+exp(-data $ x.sel %*% b.sel))
  #prob <- prob[data $ y.sel==1]
  denom <- (1+exp(-data $ x.sel %*% b.sel))
  prob <- as.matrix(1/denom[data $ y.sel==1])
  prob[is.na(prob)] <- 1e-200 # to avoid log(0), we use 1e-200
  ##===============##
  ## Weibull part
  ##===============##
  lambda <- exp( -data $ x.dur %*% b.dur)
  phi1 <- p * lambda
  phi2 <- data $ y.dur * lambda
  ft <- phi1 * (phi2^(p-1))/exp(phi2^p) ## PDF
  St <- exp(-(phi2^p))	## Survival function
  Ft <- 1-St						## CDF (1-Survival)
  ##===================================##
  ## Gaussian copula
  ##===================================##
  bound.it <- function(x,min,max) pmin(pmax(x, min), max)
  x <- bound.it(qnorm(prob),-100,100)
  y <- bound.it(qnorm(Ft),-100,100)
  J.St <- pmax(prob[t==0] - pbinorm(x[t==0], y[t==0], cov12 = b.rho), 1e-50)
  J.ft <- pmax(pnorm(q = x[t==1], mean= y[t==1]*b.rho, sd = sqrt(1-b.rho^2)) * ft[t==1], 1e-50)
  loglik <- sum(log(J.ft), log(J.St))
  return(loglik)
}

## Weibull duration with conditional logit selection

clog.weib.ll <- function(pars, data, const.vec = 0){
	## Bivariate c-logit-Weibull estimator for FIML estimation.
	## It returns log-likelihood values given data and parameters.
	##	pars = vector of parameters
	##	data = list of 4 vectors and 2 matrices.
	##   (1) selection var (y.sel)        (2) duration var (y.dur),
	##   (3) selection covariates (x.sel) (4) duration covariates (x.dur)
	##   (5) grouping var (caseid)        (6) right censoring indivator (t=0 means censoring)
	##==========================##
	## First, extract parameters
	##==========================##
	k.sel <- ncol(data $ x.sel) # n(ind v) for selection
	k.dur <- ncol(data $ x.dur) # n(ind v) for outcome
	b.sel <- pars[1:k.sel]
	b.dur <- pars[(k.sel+1):(k.sel+k.dur)]
	p <- exp(pars[(k.sel+k.dur+1)]) # log(p) -> p
	b.rho <- (exp(2*pars[length(pars)])-1)/(exp(2*pars[length(pars)])+1) # athrho -> rho
	b.rho <- min(b.rho, .9999, na.rm=T) # This is to avoid denom=0
	t <- data$t	# censoring indicator
	##===============##
	## C-logit part
	##===============##
	agg <- exp(data$x.sel %*% b.sel + const.vec) ## add constraints to the aggregator
	fac <- factor(data $ caseid)
	denom <- as.matrix(tapply(agg, fac, sum))
	prob <- agg[data $ y.sel==1]/denom
	prob[is.na(prob)] <- 1e-200 # to avoid log(0), we use 1e-200
	##===============##
	## Weibull part
	##===============##
	lambda <- exp( -data $ x.dur %*% b.dur)
	phi1 <- p * lambda
	phi2 <- data $ y.dur * lambda
	ft <- phi1 * (phi2^(p-1))/exp(phi2^p) ## PDF
	St <- exp(-(phi2^p))	## Survival function
	Ft <- 1-St						## CDF (1-Survival)
	##===================================##
	## Gaussian copula
	##===================================##
	bound.it <- function(x,min,max) pmin(pmax(x, min), max)
	x <- bound.it(qnorm(prob),-100,100)
	y <- bound.it(qnorm(Ft),-100,100)
	J.St <- pmax(prob[t==0] - pbinorm(x[t==0], y[t==0], cov12 = b.rho), 1e-50)
	J.ft <- pmax(pnorm(q = x[t==1], mean= y[t==1]*b.rho, sd = sqrt(1-b.rho^2)) * ft[t==1], 1e-50)
	loglik <- sum(log(J.ft), log(J.St))
	return(loglik)
}

## Weibull duration with logit selection (for data with no censored observation)
## Note: The only difference is that the "J.St" is 1, as pbinorm() causes
## an error when the first two arguments are numeric(0). 
## We use this function in the Monte Carlo simulation. 

log.weib.cens0.ll <- function(pars, data){
  ## Bivariate c-logit-Weibull estimator for FIML estimation.
  ## It returns log-likelihood values given data and parameters.
  ##	pars = vector of parameters
  ##	data = list of 4 vectors and 2 matrices.
  ##   (1) selection var (y.sel)        (2) duration var (y.dur),
  ##   (3) selection covariates (x.sel) (4) duration covariates (x.dur)
  ##   (5) grouping var (caseid)        (6) right censoring indivator (t=0 means censoring)
  ##==========================##
  ## First, extract parameters
  ##==========================##
  k.sel <- ncol(data $ x.sel) # n(ind v) for selection
  k.dur <- ncol(data $ x.dur) # n(ind v) for outcome
  b.sel <- pars[1:k.sel]
  b.dur <- pars[(k.sel+1):(k.sel+k.dur)]
  p <- exp(pars[(k.sel+k.dur+1)]) # log(p) -> p
  b.rho <- (exp(2*pars[length(pars)])-1)/(exp(2*pars[length(pars)])+1) # athrho -> rho
  b.rho <- min(b.rho, .9999, na.rm=T) # This is to avoid denom=0
  t <- data$t	# censoring indicator
  ##===============##
  ## logit part
  ##===============##
  denom <- (1+exp(-data $ x.sel %*% b.sel))
  prob <- 1/denom[data $ y.sel==1]
  prob[is.na(prob)] <- 1e-200 # to avoid log(0), we use 1e-200
  ##===============##
  ## Weibull part
  ##===============##
  lambda <- exp( -data $ x.dur %*% b.dur)
  phi1 <- p * lambda
  phi2 <- data $ y.dur * lambda
  ft <- phi1 * (phi2^(p-1))/exp(phi2^p) ## PDF
  St <- exp(-(phi2^p))	## Survival function
  Ft <- 1-St						## CDF (1-Survival)
  ##===================================##
  ## Gaussian copula
  ##===================================##
  bound.it <- function(x,min,max) pmin(pmax(x, min), max)
  x <- bound.it(qnorm(prob),-100,100)
  y <- bound.it(qnorm(Ft),-100,100)
  J.St <- 1
  J.ft <- pmax(pnorm(q = x[t==1], mean= y[t==1]*b.rho, sd = sqrt(1-b.rho^2)) * ft[t==1], 1e-50)
  loglik <- sum(log(J.ft), log(J.St))
  return(loglik)
}

## Weibull duration with conditional logit selection (for data with no censored observation)
## Note: The only difference is that the "J.St" is 1, as pbinorm() causes
## an error when the first two arguments are numeric(0). 
## We use this function in the Monte Carlo simulation. 

clog.weib.cens0.ll <- function(pars, data, const.vec = 0){
  ## Bivariate c-logit-Weibull estimator for FIML estimation.
  ## It returns log-likelihood values given data and parameters.
  ##	pars = vector of parameters
  ##	data = list of 4 vectors and 2 matrices.
  ##   (1) selection var (y.sel)        (2) duration var (y.dur),
  ##   (3) selection covariates (x.sel) (4) duration covariates (x.dur)
  ##   (5) grouping var (caseid)        (6) right censoring indivator (t=0 means censoring)
  ##==========================##
  ## First, extract parameters
  ##==========================##
  k.sel <- ncol(data $ x.sel) # n(ind v) for selection
  k.dur <- ncol(data $ x.dur) # n(ind v) for outcome
  b.sel <- pars[1:k.sel]
  b.dur <- pars[(k.sel+1):(k.sel+k.dur)]
  p <- exp(pars[(k.sel+k.dur+1)]) # log(p) -> p
  b.rho <- (exp(2*pars[length(pars)])-1)/(exp(2*pars[length(pars)])+1) # athrho -> rho
  b.rho <- min(b.rho, .9999, na.rm=T) # This is to avoid denom=0
  t <- data$t	# censoring indicator
  ##===============##
  ## C-logit part
  ##===============##
  agg <- exp(data$x.sel %*% b.sel + const.vec) ## add constraints to the aggregator
  fac <- factor(data $ caseid)
  denom <- as.matrix(tapply(agg, fac, sum))
  prob <- agg[data $ y.sel==1]/denom
  prob[is.na(prob)] <- 1e-200 # to avoid log(0), we use 1e-200
  ##===============##
  ## Weibull part
  ##===============##
  lambda <- exp( -data $ x.dur %*% b.dur)
  phi1 <- p * lambda
  phi2 <- data $ y.dur * lambda
  ft <- phi1 * (phi2^(p-1))/exp(phi2^p) ## PDF
  St <- exp(-(phi2^p))	## Survival function
  Ft <- 1-St						## CDF (1-Survival)
  ##===================================##
  ## Gaussian copula
  ##===================================##
  bound.it <- function(x,min,max) pmin(pmax(x, min), max)
  x <- bound.it(qnorm(prob),-100,100)
  y <- bound.it(qnorm(Ft),-100,100)
#   J.St <- pmax(prob[t==0] - pbinorm(x1 = x[t==0], x2 = y[t==0], cov12 = b.rho), 1e-50)
  J.St <- 1
  J.ft <- pmax(pnorm(q = x[t==1], mean= y[t==1]*b.rho, sd = sqrt(1-b.rho^2)) * ft[t==1], 1e-50)
  loglik <- sum(log(J.ft), log(J.St))
  return(loglik)
}

## Weibull duration with logit selection (independence assumed)

log.weib.ind.ll <- function(pars, data){
  ## Bivariate c-logit-Weibull estimator for FIML estimation (independence assumed).
  ## It returns log-likelihood values given data and parameters.
  ##	pars = vector of parameters
  ##	data = list of 4 vectors and 2 matrices.
  ##   (1) selection var (y.sel)        (2) duration var (y.dur),
  ##   (3) selection covariates (x.sel) (4) duration covariates (x.dur)
  ##   (5) grouping var (caseid)        (6) right censoring indivator (t=0 means censoring)
  ##==========================##
  ## First, extract parameters
  ##==========================##
  k.sel <- ncol(data$x.sel) # n(ind v) for selection
  k.dur <- ncol(data$x.dur) # n(ind v) for outcome
  b.sel <- pars[1:k.sel]
  b.dur <- pars[(k.sel+1):(k.sel+k.dur)]
  p <- exp(pars[(k.sel+k.dur+1)]) # log(p) -> p
  b.rho <- 0
  t <- data$t	# censoring indicator
  ##===============##
  ## logit part
  ##===============##
  denom <- (1+exp(-data $ x.sel %*% b.sel))
  prob <- 1/denom[data $ y.sel==1]
  prob[is.na(prob)] <- 1e-200 # to avoid log(0), we use 1e-200
  ##===============##
  ## Weibull part
  ##===============##
  lambda <- exp( -data$x.dur %*% b.dur)
  phi1 <- p * lambda
  phi2 <- data$y.dur * lambda
  ft <- phi1 * (phi2^(p-1))/exp(phi2^p) ## PDF
  St <- exp(-(phi2^p))	## Survival function
  Ft <- 1-St						## CDF (1-Survival)
  ##===================================##
  ## Gaussian copula
  ##===================================##
  bound.it <- function(x,min,max) pmin(pmax(x, min), max)
  x <- bound.it(qnorm(prob),-100,100)
  y <- bound.it(qnorm(Ft),-100,100)
  J.St <- pmax(prob[t==0] - pbinorm(x[t==0], y[t==0], cov12 = b.rho), 1e-50)
  J.ft <- pmax(pnorm(q = x[t==1], mean= y[t==1]*b.rho, sd = sqrt(1-b.rho^2)) * ft[t==1], 1e-50)
  loglik <- sum(log(J.ft), log(J.St))
  return(loglik)
}

## Weibull duration with conditional logit selection (independence assumed)

clog.weib.ind.ll <- function(pars, data, const.vec = 0){
	## Bivariate c-logit-Weibull estimator for FIML estimation (independence assumed).
	## It returns log-likelihood values given data and parameters.
	##	pars = vector of parameters
	##	data = list of 4 vectors and 2 matrices.
	##   (1) selection var (y.sel)        (2) duration var (y.dur),
	##   (3) selection covariates (x.sel) (4) duration covariates (x.dur)
	##   (5) grouping var (caseid)        (6) right censoring indivator (t=0 means censoring)
	##==========================##
	## First, extract parameters
	##==========================##
	k.sel <- ncol(data$x.sel) # n(ind v) for selection
	k.dur <- ncol(data$x.dur) # n(ind v) for outcome
	b.sel <- pars[1:k.sel]
	b.dur <- pars[(k.sel+1):(k.sel+k.dur)]
	p <- exp(pars[(k.sel+k.dur+1)]) # log(p) -> p
	b.rho <- 0
	t <- data$t	# censoring indicator
	##===============##
	## C-logit part
	##===============##
	agg <- exp(data$x.sel %*% b.sel + const.vec) ## add constraints to the aggregator
	fac <- factor(data$caseid)
	denom <- as.matrix(tapply(agg, fac, sum))
	prob <- agg[data$y.sel==1]/denom
	prob[is.na(prob)] <- 1e-200 # to avoid log(0), we use 1e-200
	##===============##
	## Weibull part
	##===============##
	lambda <- exp( -data$x.dur %*% b.dur)
	phi1 <- p * lambda
	phi2 <- data$y.dur * lambda
	ft <- phi1 * (phi2^(p-1))/exp(phi2^p) ## PDF
	St <- exp(-(phi2^p))	## Survival function
	Ft <- 1-St						## CDF (1-Survival)
	##===================================##
	## Gaussian copula
	##===================================##
	bound.it <- function(x,min,max) pmin(pmax(x, min), max)
	x <- bound.it(qnorm(prob),-100,100)
	y <- bound.it(qnorm(Ft),-100,100)
	J.St <- pmax(prob[t==0] - pbinorm(x1 = x[t==0], x2 = y[t==0], cov12 = b.rho), 1e-50)
	J.ft <- pmax(pnorm(q = x[t==1], mean= y[t==1]*b.rho, sd = sqrt(1-b.rho^2)) * ft[t==1], 1e-50)
	loglik <- sum(log(J.ft), log(J.St))
	return(loglik)
}

## Log-logistic duration with conditional logit selection

clog.loglog.ll <- function(pars, data, const.vec = 0){
	## Bivariate c-logit-loglogistic estimator for FIML estimation.
	## It returns log-likelihood values given data and parameters.
	##	pars = vector of parameters
	##	data = list of 4 vectors and 2 matrices.
	##   (1) selection var (y.sel)        (2) duration var (y.dur),
	##   (3) selection covariates (x.sel) (4) duration covariates (x.dur)
	##   (5) grouping var (caseid)        (6) right censoring indivator (t=0 means censoring)
	##==========================##
	## First, extract parameters
	##==========================##
	k.sel <- ncol(data$x.sel) # n(ind v) for selection
	k.dur <- ncol(data$x.dur) # n(ind v) for outcome
	b.sel <- pars[1:k.sel]
	b.dur <- pars[(k.sel+1):(k.sel+k.dur)]
	p <- exp(pars[(k.sel+k.dur+1)]) # log(p) -> p
	b.rho <- (exp(2*pars[length(pars)])-1)/(exp(2*pars[length(pars)])+1) # athrho -> rho
	b.rho <- min(b.rho, .9999, na.rm=T) # This is to avoid denom=0
	t <- data$t	# censoring indicator
	##===============##
	## C-logit part
	##===============##
	agg <- exp(data$x.sel %*% b.sel + const.vec) ## add constraints to the aggregator
	fac <- factor(data$caseid)
	denom <- as.matrix(tapply(agg, fac, sum))
	prob <- agg[data$y.sel==1]/denom
	prob[is.na(prob)] <- 1e-200 # to avoid log(0), we use 1e-200
	##===============##
	## Loglogit part
	##===============##
	lambda <- exp( -data$x.dur %*% b.dur )
	ft <- ( lambda * p * (lambda * data$y.dur)^(p-1) )/(( 1 + ( lambda * data$y.dur )^p )^2)	
	St <- 1/(1 + ( lambda * data$y.dur )^p )
	Ft <- 1-St						## CDF (1-Survival)
	##===================================##
	## Gaussian copula
	##===================================##
	bound.it <- function(x,min,max) pmin(pmax(x, min), max)
	x <- bound.it(qnorm(prob),-100,100)
	y <- bound.it(qnorm(Ft),-100,100)
	J.St <- pmax(prob[t==0] - pbinorm(x[t==0], y[t==0], cov12 = b.rho), 1e-50)
	J.ft <- pmax(pnorm(q = x[t==1], mean= y[t==1]*b.rho, sd = sqrt(1-b.rho^2)) * ft[t==1], 1e-50)
	loglik <- sum(log(J.ft), log(J.St))
	return(loglik)
}

## Log-logistic duration with conditional logit selection (independence assumed)

clog.loglog.ind.ll <- function(pars, data, const.vec = 0){
	## Bivariate c-logit-Loglogistic estimator for FIML estimation (independence assumed).
	## It returns log-likelihood values given data and parameters.
	##	pars = vector of parameters
	##	data = list of 4 vectors and 2 matrices.
	##   (1) selection var (y.sel)        (2) duration var (y.dur),
	##   (3) selection covariates (x.sel) (4) duration covariates (x.dur)
	##   (5) grouping var (caseid)        (6) right censoring indivator (t=0 means censoring)
	##==========================##
	## First, extract parameters
	##==========================##
	k.sel <- ncol(data$x.sel) # n(ind v) for selection
	k.dur <- ncol(data$x.dur) # n(ind v) for outcome
	b.sel <- pars[1:k.sel]
	b.dur <- pars[(k.sel+1):(k.sel+k.dur)]
	p <- exp(pars[(k.sel+k.dur+1)]) # log(p) -> p
	b.rho <- 0
	t <- data$t	# censoring indicator
	##===============##
	## C-logit part
	##===============##
	agg <- exp(data$x.sel %*% b.sel + const.vec) ## add constraints to the aggregator
	fac <- factor(data$caseid)
	denom <- as.matrix(tapply(agg, fac, sum))
	prob <- agg[data$y.sel==1]/denom
	prob[is.na(prob)] <- 1e-200 # to avoid log(0), we use 1e-200
	##===============##
	## Weibull part
	##===============##
	lambda <- exp( -data$x.dur %*% b.dur)
	phi1 <- p * lambda
	phi2 <- data$y.dur * lambda
	ft <- ( lambda * p * (lambda * data$y.dur)^(p-1) )/(( 1 + ( lambda * data$y.dur )^p )^2)	
	St <- 1/(1 + ( lambda * data$y.dur )^p )
	Ft <- 1-St						## CDF (1-Survival)
	##===================================##
	## Gaussian copula
	##===================================##
	bound.it <- function(x,min,max) pmin(pmax(x, min), max)
	x <- bound.it(qnorm(prob),-100,100)
	y <- bound.it(qnorm(Ft),-100,100)
	J.St <- pmax(prob[t==0] - pbinorm(x1 = x[t==0], x2 = y[t==0], cov12 = b.rho), 1e-50)
	J.ft <- pmax(pnorm(q = x[t==1], mean= y[t==1]*b.rho, sd = sqrt(1-b.rho^2)) * ft[t==1], 1e-50)
	loglik <- sum(log(J.ft), log(J.St))
	return(loglik)
}


# Auxiliary functions for ML implementation -------------------------------

## function to convert ath(rho) to rho

get.rho <- function(x) (exp(2*x)-1)/(exp(2*x)+1)

## wrapper function for log.weib.ll

estimate.ld <- function(
  eqn.sel, eqn.dur,	## selection & duration equations  (formula class)
  failure.type, 		## failure/censoring indicator (=1 if failed, =0 if censored)
  dist = "",        ## distribution ("weibull" or "loglogistic")
  data,             ## Data (data frame class)
  WantHessian = T,  ## calculate Hessian as well? (logical)
  init  				    ## initial values
){
  
  if (is.null(data)) stop("No data provided")
  if (dist=="weibull") {dist.ll <- log.weib.ll} else{
    if (dist=="loglogistic") {dist.ll <- clog.loglog.ll} else{
      stop("Must specify an appropriate distribution")}}
  if (is.null(data)) stop("No data provided")
  if (failure.type=="pooled") {failure <- data$censor_pooled} else{
    if (failure.type=="replacement") {failure <- data$censor_replace} else{
      if (failure.type=="election") {failure <- data$censor_discelec} else{
        stop("Must specify an appropriate failure type")}}}

  # extract variables
  name.of.y.sel <- rownames(attr(terms(eqn.sel),"factors"))[1]
  name.of.y.dur <- rownames(attr(terms(eqn.dur),"factors"))[1]
  names.of.x.sel <- attr(terms(eqn.sel),"term.labels")
  names.of.x.dur <- attr(terms(eqn.dur),"term.labels")
  y.sel <- data[, match(name.of.y.sel, colnames(data))]
  y.dur <- data[y.sel==1, match(name.of.y.dur, colnames(data))]
  x.sel <- as.matrix(data[, match(names.of.x.sel,colnames(data))])
  x.sel <- cbind(x.sel,1)
  x.dur <- as.matrix(data[y.sel==1, match(names.of.x.dur,colnames(data))])
  x.dur <- cbind(x.dur,1) ## add intercept
  t <- failure[y.sel==1]
  t[is.na(y.dur)==1] <- 2	
  estim.data <- list(y.sel = y.sel, y.dur = y.dur, x.sel = x.sel, x.dur = x.dur, t = t)
  # Fit the model
  fit <- optim(par = init, ## initial values
               fn = dist.ll,
               data = estim.data,
               method = "BFGS", ## BFGS, Nelder-Mead, SANN, CG
               control = list(fnscale = -1, maxit = 10000, trace=1), 
               hessian = WantHessian)
  # Return
  out <- list(fit = fit, sel.covs = names.of.x.sel, dur.covs = names.of.x.dur)
  return(out)
}

## wrapper function for clog.weib.ll and clog.loglog.ll

estimate.cld <- function(
    eqn.sel, eqn.dur,	## selection & duration equations  (formula class)
    failure.type, 		## failure/censoring indicator (=1 if failed, =0 if censored)
    dist = "",        ## distribution ("weibull" or "loglogistic")
    data,             ## Data (data frame class)
    WantHessian = T,  ## calculate Hessian as well? (logical)
    init, 				    ## initial values
    constraint = T    ## use constraint? (logical)
){
  
  if (is.null(data)) stop("No data provided")
  if (dist=="weibull") {dist.ll <- clog.weib.ll} else{
    if (dist=="loglogistic") {dist.ll <- clog.loglog.ll} else{
      stop("Must specify an appropriate distribution")}}
  if (is.null(data)) stop("No data provided")
  if (failure.type=="pooled") {failure <- data$censor_pooled} else{
    if (failure.type=="replacement") {failure <- data$censor_replace} else{
      if (failure.type=="election") {failure <- data$censor_discelec} else{
        stop("Must specify an appropriate failure type")}}}
  
  # constraint
  if (constraint == T){
    const.vec.n <- as.matrix(data[const.neg]) %*% rep(-3, length(const.neg))
    const.vec.p <- as.matrix(data[const.pos]) %*% rep( 3, length(const.pos))
    const.vec <- const.vec.n + const.vec.p    
  }
  if (constraint == F) const.vec <- 0
  
  # extract variables
  name.of.y.sel <- rownames(attr(terms(eqn.sel),"factors"))[1]
  name.of.y.dur <- rownames(attr(terms(eqn.dur),"factors"))[1]
  names.of.x.sel <- attr(terms(eqn.sel),"term.labels")
  names.of.x.dur <- attr(terms(eqn.dur),"term.labels")
  y.sel <- data[, match(name.of.y.sel, colnames(data))]
  y.dur <- data[y.sel==1, match(name.of.y.dur, colnames(data))]
  x.sel <- as.matrix(data[, match(names.of.x.sel,colnames(data))])
  x.dur <- as.matrix(data[y.sel==1, match(names.of.x.dur,colnames(data))])
  x.dur <- cbind(x.dur,1) ## add intercept
  t <- failure[y.sel==1]
  t[is.na(y.dur)==1] <- 2	
  estim.data <- list(y.sel = y.sel, y.dur = y.dur, x.sel = x.sel, x.dur = x.dur,
                     caseid = data$caseid, t = t)
  # Fit the model
  fit <- optim(par = init, ## initial values
               fn = dist.ll,
               data = estim.data,
               const.vec = const.vec, 
               method = "SANN", ## BFGS, Nelder-Mead, SANN, CG
               control = list(fnscale = -1, maxit = 30000, trace=1), 
               hessian = WantHessian)
  # Return
  out <- list(fit = fit, sel.covs = names.of.x.sel, dur.covs = names.of.x.dur)
  return(out)
}

## Function to display the output from estimate.cld function

disp.cld <- function(estim, rnd=4){
	fit <- estim$fit
	beta <- fit$par
	vcov <- -solve(fit$hessian)
	vcov <- make.positive.definite(vcov)
	se <- sqrt(diag(vcov))
	## transform ath(rho) to rho
	beta <- c(beta, get.rho(beta[length(beta)]))
	se <- c(se, NA)
	## transform ln(p) to p
	beta <- c(beta, exp(beta[(length(beta)-2)]))
	se <- c(se, NA)
	## z-score, p-value
	z <- beta/se
	pval <- 2*(1-pnorm(abs(beta/se)))
	tbl <- cbind(beta = round(beta, rnd),
							 se = round(se, rnd),
							 z = round(z,3),
							 pval = round(pval, 3))
	## separate by stages
	k.sel <- length(estim$sel.covs)
	k.dur <- length(estim$dur.covs) + 1
	b.sel <- tbl[1:k.sel,]									 # Beta for selection
	b.dur <- tbl[(k.sel+1):(k.sel+k.dur+1),] # Beta for duration
	b.dur <- rbind(b.dur, tbl[(k.sel+k.dur+4),]) # add exp(ln(p))=p
	b.rho <- tbl[(k.sel+k.dur+2):(k.sel+k.dur+3),] # Rho
	## var names
	rownames(b.sel) <- estim$sel.covs
	rownames(b.dur) <- c(estim$dur.covs, "Constant","log(p)","p")
	rownames(b.rho) <- c("ath(rho)","rho")
	## bind
	out <- list(Selection = b.sel,
							Duration = b.dur,
							Correlation = b.rho)
	return(out)
}

## wrapper function for log.weib.ind.ll

estimate.ild <- function(
  eqn.sel, eqn.dur,  ## selection & duration equations  (formula class)
  failure.type, 		## failure/censoring indicator (=1 if failed, =0 if censored)
  dist = "",        ## distribution ("weibull" or "loglogistic")
  data,             ## Data (data frame class)
  WantHessian = T,  ## calculate Hessian as well? (logical)
  init 				    ## initial values
){
  
  if (is.null(data)) stop("No data provided")
  if (dist=="weibull") {dist.ll <- log.weib.ind.ll} else{
    if (dist=="loglogistic") {dist.ll <- clog.loglog.ind.ll} else{
      stop("Must specify an appropriate distribution")}}
  if (is.null(data)) stop("No data provided")
  if (failure.type=="pooled") {failure <- data$censor_pooled} else{
    if (failure.type=="replacement") {failure <- data$censor_replace} else{
      if (failure.type=="election") {failure <- data$censor_discelec} else{
        stop("Must specify an appropriate failure type")}}}
  
  # extract variables
  name.of.y.sel <- rownames(attr(terms(eqn.sel),"factors"))[1]
  name.of.y.dur <- rownames(attr(terms(eqn.dur),"factors"))[1]
  names.of.x.sel <- attr(terms(eqn.sel),"term.labels")
  names.of.x.dur <- attr(terms(eqn.dur),"term.labels")
  y.sel <- data[, match(name.of.y.sel, colnames(data))]
  y.dur <- data[y.sel==1, match(name.of.y.dur, colnames(data))]
  x.sel <- as.matrix(data[, match(names.of.x.sel,colnames(data))])
  x.dur <- as.matrix(data[y.sel==1, match(names.of.x.dur,colnames(data))])
  x.dur <- cbind(x.dur,1) ## add intercept
  t <- failure[y.sel==1]
  t[is.na(y.dur)==1] <- 2
  estim.data <- list(y.sel = y.sel, y.dur = y.dur, x.sel = x.sel, x.dur = x.dur,
                     caseid = data$caseid, t = t)
  # Fit the model
  fit <- optim(par = init, ## initial values
               fn = dist.ll,
               data = estim.data,
               method = "BFGS", ## BFGS, Nelder-Mead, SANN, CG
               control = list(fnscale = -1, maxit = 10000, trace=1),
               hessian = WantHessian)
  # Return
  out <- list(fit = fit, sel.covs = names.of.x.sel, dur.covs = names.of.x.dur)
  return(out)
}

## wrapper function for log.weib.cens0.ll

estimate.ldcens0 <- function(
  eqn.sel, eqn.dur,	## selection & duration equations  (formula class)
  failure.type, 		## failure/censoring indicator (=1 if failed, =0 if censored)
  dist = "",        ## distribution ("weibull" or "loglogistic")
  data,             ## Data (data frame class)
  WantHessian = T,  ## calculate Hessian as well? (logical)
  init  				    ## initial values
){
  
  if (is.null(data)) stop("No data provided")
  if (dist=="weibull") {dist.ll <- log.weib.cens0.ll} else{
    if (dist=="loglogistic") {dist.ll <- clog.loglog.ll} else{
      stop("Must specify an appropriate distribution")}}
  if (is.null(data)) stop("No data provided")
  if (failure.type=="pooled") {failure <- data$censor_pooled} else{
    if (failure.type=="replacement") {failure <- data$censor_replace} else{
      if (failure.type=="election") {failure <- data$censor_discelec} else{
        stop("Must specify an appropriate failure type")}}}
  
  # extract variables
  name.of.y.sel <- rownames(attr(terms(eqn.sel),"factors"))[1]
  name.of.y.dur <- rownames(attr(terms(eqn.dur),"factors"))[1]
  names.of.x.sel <- attr(terms(eqn.sel),"term.labels")
  names.of.x.dur <- attr(terms(eqn.dur),"term.labels")
  y.sel <- data[, match(name.of.y.sel, colnames(data))]
  y.dur <- data[y.sel==1, match(name.of.y.dur, colnames(data))]
  x.sel <- as.matrix(data[, match(names.of.x.sel,colnames(data))])
  x.sel <- cbind(x.sel,1)
  x.dur <- as.matrix(data[y.sel==1, match(names.of.x.dur,colnames(data))])
  x.dur <- cbind(x.dur,1) ## add intercept
  t <- failure[y.sel==1]
  t[is.na(y.dur)==1] <- 2	
  estim.data <- list(y.sel = y.sel, y.dur = y.dur, x.sel = x.sel, x.dur = x.dur, t = t)
  # Fit the model
  fit <- optim(par = init, ## initial values
               fn = dist.ll,
               data = estim.data,
               method = "BFGS", ## BFGS, Nelder-Mead, SANN, CG
               control = list(fnscale = -1, maxit = 10000, trace=1), 
               hessian = WantHessian)
  # Return
  out <- list(fit = fit, sel.covs = names.of.x.sel, dur.covs = names.of.x.dur)
  return(out)
}

## wrapper function for clog.weib.ind.ll and clog.loglog.ind.ll

estimate.icld <- function(
  eqn.sel, eqn.dur,  ## selection & duration equations  (formula class)
  failure.type, 		## failure/censoring indicator (=1 if failed, =0 if censored)
  dist = "",        ## distribution ("weibull" or "loglogistic")
  data,             ## Data (data frame class)
  WantHessian = T,  ## calculate Hessian as well? (logical)
  init, 				    ## initial values
  constraint = T    ## use constraint? (logical)
){
  
  if (is.null(data)) stop("No data provided")
  if (dist=="weibull") {dist.ll <- clog.weib.ind.ll} else{
    if (dist=="loglogistic") {dist.ll <- clog.loglog.ind.ll} else{
      stop("Must specify an appropriate distribution")}}
  if (is.null(data)) stop("No data provided")
  if (failure.type=="pooled") {failure <- data$censor_pooled} else{
    if (failure.type=="replacement") {failure <- data$censor_replace} else{
      if (failure.type=="election") {failure <- data$censor_discelec} else{
        stop("Must specify an appropriate failure type")}}}
  
  # constraint
  if (constraint == T){
    const.vec.n <- as.matrix(data[const.neg]) %*% rep(-3, length(const.neg))
    const.vec.p <- as.matrix(data[const.pos]) %*% rep( 3, length(const.pos))
    const.vec <- const.vec.n + const.vec.p    
  }
  if (constraint == F) const.vec <- 0
  
  # extract variables
  name.of.y.sel <- rownames(attr(terms(eqn.sel),"factors"))[1]
  name.of.y.dur <- rownames(attr(terms(eqn.dur),"factors"))[1]
  names.of.x.sel <- attr(terms(eqn.sel),"term.labels")
  names.of.x.dur <- attr(terms(eqn.dur),"term.labels")
  y.sel <- data[, match(name.of.y.sel, colnames(data))]
  y.dur <- data[y.sel==1, match(name.of.y.dur, colnames(data))]
  x.sel <- as.matrix(data[, match(names.of.x.sel,colnames(data))])
  x.dur <- as.matrix(data[y.sel==1, match(names.of.x.dur,colnames(data))])
  x.dur <- cbind(x.dur,1) ## add intercept
  t <- failure[y.sel==1]
  t[is.na(y.dur)==1] <- 2
  estim.data <- list(y.sel = y.sel, y.dur = y.dur, x.sel = x.sel, x.dur = x.dur,
                     caseid = data$caseid, t = t)
	# Fit the model
	fit <- optim(par = init, ## initial values
							 fn = dist.ll,
							 data = estim.data,
							 const.vec = const.vec, 
							 method = "BFGS", ## BFGS, Nelder-Mead, SANN, CG
							 control = list(fnscale = -1, maxit = 10000, trace=1),
							 hessian = WantHessian)
	# Return
	out <- list(fit = fit, sel.covs = names.of.x.sel, dur.covs = names.of.x.dur)
	return(out)
}

## Function to display the output from estimate.icld function

disp.icld <- function(estim, rnd=4){
	fit <- estim$fit
	beta <- fit$par
	vcov <- -solve(fit$hessian)
	vcov <- make.positive.definite(vcov)
	se <- sqrt(diag(vcov))
	## transform ln(p) to p
	beta <- c(beta, exp(beta[(length(beta))]))
	se <- c(se, NA)
	## z-score, p-value
	z <- beta/se
	pval <- 2*(1-pnorm(abs(beta/se)))
	tbl <- cbind(beta = round(beta, rnd),
							 se = round(se, rnd),
							 z = round(z,3),
							 pval = round(pval, 3))
	## separate by stages
	k.sel <- length(estim$sel.covs)
	k.dur <- length(estim$dur.covs) + 1
	b.sel <- tbl[1:k.sel,]									 # Beta for selection
	b.dur <- tbl[(k.sel+1):(k.sel+k.dur+1),] # Beta for duration
	b.dur <- rbind(b.dur, tbl[(k.sel+k.dur+2),]) # add exp(ln(p))=p
	## var names
	rownames(b.sel) <- estim$sel.covs
	rownames(b.dur) <- c(estim$dur.covs, "Constant","log(p)","p")
	## bind
	out <- list(Selection = b.sel,
							Duration = b.dur,
							Correlation = 0)
	return(out)
}

## wrapper function for logit.ll

estimate.logit.prob <- function(
  eqn.sel,          ## Selection equation (formula class)
  data,             ## Data (data frame class)
  WantHessian = T  ## calculate Hessian as well? (logical)
){
  
  # extract variables
  name.of.y.sel <- rownames(attr(terms(eqn.sel),"factors"))[1]
  names.of.x.sel <- attr(terms(eqn.sel),"term.labels")
  y.sel <- as.matrix(data[, match(name.of.y.sel, colnames(data))])
  x.sel <- as.matrix(data[, match(names.of.x.sel,colnames(data))])
  x.sel <- cbind(x.sel,1)
  estim.data <- list(y.vec = y.sel, x.mat = x.sel)
  
  # Fit the model
  par.init.val = rep(0,ncol(x.sel)) # arbitrary starting parameters
  fit <- optim(par = par.init.val, ## initial values
               fn = logit.ll.prob,
               data = estim.data,
               method = "BFGS", 
               hessian = WantHessian)
  # Return
  out <- list(fit = fit, sel.covs = names.of.x.sel)
  return(out)
}

## wrapper function for clogit.ll

estimate.clogit <- function(
  eqn.sel,          ## Selection equation (formula class)
  data,             ## Data (data frame class)
  WantHessian = T,  ## calculate Hessian as well? (logical)
  constraint = T    ## use constraint? (logical)
  ){
  # constraint
  if (constraint == T){
    const.vec.n <- as.matrix(data[const.neg]) %*% rep(-3, length(const.neg))
    const.vec.p <- as.matrix(data[const.pos]) %*% rep( 3, length(const.pos))
    const.vec <- const.vec.n + const.vec.p    
  }
  if (constraint == F) const.vec <- 0
  
  # extract variables
	name.of.y.sel <- rownames(attr(terms(eqn.sel),"factors"))[1]
	names.of.x.sel <- attr(terms(eqn.sel),"term.labels")
	y.sel <- data[, match(name.of.y.sel, colnames(data))]
	x.sel <- as.matrix(data[, match(names.of.x.sel,colnames(data))])
	estim.data <- list(y = y.sel, xmat = x.sel, caseid = data$caseid)

  # Fit the model
	k <- ncol(x.sel)
	fit <- optim(par = rep(0,k), ## initial values
							 fn = clogit.ll,
							 data = estim.data,
               const.vec = const.vec,
							 method = "BFGS", 
							 control = list(fnscale = -1, maxit = 10000, trace=1),
							 hessian = WantHessian)
	# Return
	out <- list(fit = fit,
							sel.covs = names.of.x.sel)
	return(out)
}

## Function to display the output from estimate.clogit function

disp.clogit <- function(estim, rnd=4){
	fit <- estim$fit
	beta <- fit$par
	vcov <- -solve(fit$hessian)
	vcov <- make.positive.definite(vcov)
	se <- sqrt(diag(vcov))
	## z-score, p-value
	z <- beta/se
	pval <- 2*(1-pnorm(abs(beta/se)))
	tbl <- cbind(beta = round(beta, rnd),
							 se = round(se, rnd),
							 z = round(z,3),
							 pval = round(pval, 3))
	## separate by stages
	k <- length(estim$sel.covs)
	## var names
	rownames(tbl) <- estim$sel.covs
	out <- list(Selection = tbl)
	return(out)
}

## wrapper function for weibull.ll and loglogistic.ll

estimate.univdur <- function(
  eqn.dur, ## Duration equation  (formula class)
  failure.type,   	## failure/censoring indicator (=1 if failed, =0 if censored)
  dist = "",        ## distribution ("weibull" or "loglogistic")
  data              ## Data (data frame class)
  ){
	if (is.null(data)) stop("No data provided")
	if (dist=="weibull") {dist.ll <- weibull.ll} else{
	if (dist=="loglogistic") {dist.ll <- loglogistic.ll} else{
	stop("Must specify appropriate distributions")}}
	if (failure.type=="pooled") {failure <- data$censor_pooled} else{
	if (failure.type=="replacement") {failure <- data$censor_replace} else{
	if (failure.type=="election") {failure <- data$censor_discelec} else{
	stop("Must specify appropriate failure type")}}}

	name.of.y.dur <- rownames(attr(terms(eqn.dur),"factors"))[1]
	names.of.x.dur <- attr(terms(eqn.dur),"term.labels")
	y.dur <- data[, match(name.of.y.dur, colnames(data))]
	id <- is.na(y.dur)!=1
	y.dur <- y.dur[id]
	x.dur <- as.matrix(data[id, match(names.of.x.dur,colnames(data))])
	x.dur <- cbind(x.dur,1) ## add intercept
	t <- failure[id]

	estim.data <- list(y = y.dur, xmat = x.dur, t = t)

  # Fit the model
	k <- ncol(x.dur) + 1 ## add ln(p)
	fit <- optim(par = rep(0.1,k), ## initial values
							 fn = dist.ll,
							 data = estim.data,
							 method = "BFGS", ## BFGS, Nelder-Mead, SANN, CG
							 control = list(fnscale = -1, maxit = 10000, trace=1),
							 hessian = T)

  # Return
	out <- list(fit = fit,
							dur.covs = names.of.x.dur)
	return(out)
}


## wrapper function for weibull.ll.str

estimate.univdur.str <- function(
  eqn.dur, ## Duration equation  (formula class)
  failure.type,   	## failure/censoring indicator (=1 if failed, =0 if censored)
  dist = "",        ## distribution ("weibull" or "loglogistic")
  data              ## Data (data frame class)
){
  if (is.null(data)) stop("No data provided")
  if (dist=="weibull") {dist.ll <- weibull.ll.str} else{
    if (dist=="loglogistic") {dist.ll <- loglogistic.ll} else{
      stop("Must specify appropriate distributions")}}
  if (failure.type=="pooled") {failure <- data$censor_pooled} else{
    if (failure.type=="replacement") {failure <- data$censor_replace} else{
      if (failure.type=="election") {failure <- data$censor_discelec} else{
        stop("Must specify appropriate failure type")}}}
  
  name.of.y.dur <- rownames(attr(terms(eqn.dur),"factors"))[1]
  names.of.x.dur <- attr(terms(eqn.dur),"term.labels")
  y.dur <- data[, match(name.of.y.dur, colnames(data))]
  id <- is.na(y.dur)!=1
  y.dur <- y.dur[id]
  x.dur <- as.matrix(data[id, match(names.of.x.dur,colnames(data))])
  x.dur <- cbind(x.dur,1) ## add intercept
  t <- failure[id]
  strata <- as.matrix(data[id, 40])
  estim.data <- list(y = y.dur, xmat = x.dur, t = t, strata = strata)
  
  # Fit the model
  k <- ncol(x.dur) + 1 ## add ln(p)
  fit <- optim(par = rep(0.1,k), ## initial values
               fn = dist.ll,
               data = estim.data,
               method = "BFGS", ## BFGS, Nelder-Mead, SANN, CG
               control = list(fnscale = -1, maxit = 10000, trace=1),
               hessian = T)
  
  # Return
  out <- list(fit = fit,
              dur.covs = names.of.x.dur)
  return(out)
}

# Function to display the output from estimate.univdur

disp.univdur <- function(estim, rnd=4){
	fit <- estim$fit
	beta <- fit$par
	vcov <- -solve(fit$hessian)
	vcov <- make.positive.definite(vcov)
	se <- sqrt(diag(vcov))
	## transform ln(p) to p
	beta <- c(beta, exp(beta[length(beta)]))
	se <- c(se, NA)
	## z-score, p-value
	z <- beta/se
	pval <- 2*(1-pnorm(abs(beta/se)))
	tbl <- cbind(beta = round(beta, rnd),
							 se = round(se, rnd),
							 z = round(z,3),
							 pval = round(pval, 3))
	## var names
	rownames(tbl) <- c(estim$dur.covs, "Constant","log(p)","p")
	## bind
	out <- list(Duration = tbl)
	return(out)
}


hmftest.clogit <- function(model, drop_proportion ,iter, ...){
  if(missing(drop_proportion)) {
    drop_proportion <- 0.1
    message("No drop proportion specified, so 0.1 was chosen")
  }
  if(missing(iter)) {
    iter <- 1
    message("No numer of iterations was specified, so 1 was chosen")
  }
  
  formulastring <- as.character(model$userCall)[2]
  formula <- as.formula(formulastring)
  dataused <- as.character(model$userCall)[3]
  dataused <- get(dataused)
  strata <- sub(".*strata\\(([^\\)]+)\\).*", "\\1", formulastring)
  depvar <- sub("^(\\S+).*", "\\1", formulastring)
  
  
  if (iter==1) {
    temp_data <- dataused %>%
      group_by(!!sym(strata)) %>%
      mutate(
        can_drop = ifelse(!!sym(depvar) == 1, FALSE, TRUE)
      ) %>%
      filter(
        !!sym(depvar) == 1 | {
          can_drop_values <- which(can_drop)  
          if (length(can_drop_values) > 0) {  
            num_to_drop <- max(1, ceiling((1-drop_proportion) * length(can_drop_values)))
            rows_to_drop <- sample(can_drop_values, size = num_to_drop, replace = FALSE)
            row_number() %in% rows_to_drop
          } else {
            FALSE 
          }
        }
      ) %>%
      ungroup() %>%
      select(-can_drop)
    
    altmodel <- clogit(formula, temp_data)
    
    coef.x <- coef(model)
    coef.s <- coef(altmodel)
    un <- names(coef.x) %in% names(coef.s)
    diff.coef <- coef.s - coef.x[un]
    diff.var <- vcov(altmodel) - vcov(model)[un, un]
    hmf <- as.numeric(diff.coef %*% solve(diff.var) %*% diff.coef)
    names(hmf) <- "chisq"
    df <- sum(un)
    names(df) <- "df"
    pv <- pchisq(hmf, df = df, lower.tail = FALSE)
    res <- list(data.name = model$call$data,
                statistic = hmf,
                p.value =pv,
                parameter = df,
                method = "Hausman-McFadden test",
                alternative = "IIA is rejected")
    class(res) <- "htest"
    res
  } else {
    output <- data.frame(pval = numeric(0), chisq = numeric(0))
    for (iterations in 1:iter) {
      temp_data <- dataused %>%
        group_by(!!sym(strata)) %>%
        mutate(
          can_drop = ifelse(!!sym(depvar) == 1, FALSE, TRUE)
        ) %>%
        filter(
          !!sym(depvar) == 1 | {
            can_drop_values <- which(can_drop)  
            if (length(can_drop_values) > 0) {  
              num_to_drop <- max(1, ceiling((1-drop_proportion) * length(can_drop_values)))
              rows_to_drop <- sample(can_drop_values, size = num_to_drop, replace = FALSE)
              row_number() %in% rows_to_drop
            } else {
              FALSE 
            }
          }
        ) %>%
        ungroup() %>%
        select(-can_drop)
      
      altmodel <- clogit(formula, temp_data)
      
      coef.x <- coef(model)
      coef.s <- coef(altmodel)
      un <- names(coef.x) %in% names(coef.s)
      diff.coef <- coef.s - coef.x[un]
      diff.var <- vcov(altmodel) - vcov(model)[un, un]
      hmf <- as.numeric(diff.coef %*% solve(diff.var) %*% diff.coef)
      df <- sum(un)
      pv <- pchisq(hmf, df = df, lower.tail = FALSE)  
      output[iterations,1] <- pv
      output[iterations,2] <- hmf
    }
    output
  }
}  







# End of file -------------------------------------------------------------
